Report no. 07/18 



Kalman Filtering with Equality and Inequality State 

Constraints 



Nachi Gupta Raphael Hauser 

Oxford University Computing Laboratory, Numerical Analysis Group, 
Wolfson Building, Parks Road, Oxford 0X1 3QD, U.K. 



Both constrained and unconstrained optimization problems regularly appear 
in recursive tracking problems engineers currently address - however, constraints 
are rarely exploited for these applications. We define the Kalman Filter and 
discuss two different approaches to incorporating constraints. Each of these ap- 
proaches are first applied to equality constraints and then extended to inequality 
constraints. We discuss methods for dealing with nonlinear constraints and for 
constraining the state prediction. Finally, some experiments are provided to indi- 
cate the usefulness of such methods. 



Key words and phrases: Constrained optimization, Kalman filtering. Nonlinear filters. 
Optimization methods. Quadratic programming. State estimation 

The first author would like to thank the Clarendon Bursary for financial support. 



Oxford University Computing Laboratory 
Numerical Analysis Group 
Wolfson Building 
Parks Road 

Oxford, England 0X1 3QD 

E-mail: nachi@comlab.oxford.ac.uk February, 2008 



2 



1 Introduction 



Kalman Filtering [8] is a method to make real-time predictions for systems with some known 
dynamics. Traditionally, problems requiring Kalman Filtering have been complex and non- 
linear. Many advances have been made in the direction of dealing with nonlinearities (e.g., 
Extended Kalman Filter [1], Unscented Kalman Filter [7]). These problems also tend to have 
inherent state space equality constraints (e.g., a fixed speed for a robotic arm) and state space 
inequality constraints (e.g., maximum attainable speed of a motor). In the past, less interest 
has been generated towards constrained Kalman Filtering, partly because constraints can be 
difficult to model. As a result, constraints are often neglected in standard Kalman Filtering 
applications. 

The extension to Kalman Filtering with known equality constraints on the state space is 
discussed in [5, 11, 13, 14, 16]. In this paper, we discuss two distinct methods to incorporate 
constraints into a Kalman Filter. Initially, we discuss these in the framework of equality con- 
straints. The first method, projecting the updated state estimate onto the constrained region, 
appears with some discussion in [5, 1 1]. We propose another method, which is to restrict the 
optimal Kalman Gain so the updated state estimate will not violate the constraint. With some 
algebraic manipulation, the second method is shown to be a special case of the first method. 

We extend both of these concepts to Kalman Filtering with inequality constraints in the 
state space. This generalization for the first approach was discussed in [12]Q Constraining the 
optimal Kalman Gain was briefly discussed in [10]. Further, we will also make the extension 
to incorporating state space constraints in Kalman Filter predictions. 

Analogous to the way a Kalman Filter can be extended to solve problems containing non- 
linearities in the dynamics using an Extended Kalman Filter by linearizing locally (or by using 
an Unscented Kalman Filter), linear inequality constrained filtering can similarly be extended 
to problems with nonlinear constraints by linearizing locally (or by way of another scheme 
like an Unscented Kalman Filter). The accuracy achieved by methods dealing with nonlinear 
constraints will naturally depend on the structure and curvature of the nonlinear function it- 
self. In the two experiments we provide, we look at incorporating inequality constraints to a 
tracking problem with nonlinear dynamics. 



2 Kalman Filter 

A discrete-time Kalman Filter [8] attempts to find the best running estimate for a recursive 
system governed by the following modeo 

Xk = Fk,k-iXk-i+Uk,k-i-: Uk,k-i ^ ^ {0,Qk,k-i) (2.1) 



Zk = HkXk + Vk, Vk^^ (0, Rk) (2.2) 
'The similar extension for the method of [16] was made in [6]. 

^The subscript A: on a variable stands for the fc-th time step, the mathematical notation .VV {jJ.,!.) denotes a 
normally distributed random vector with mean jj, and covariance E, and all vectors in this paper are column 
vectors (unless we are explicitly taking the transpose of the vector). 
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Here xi^ is an n-vector that represents the true state of the underlying system and Fj^ j^^i 
is an n X n matrix that describes the transition dynamics of the system from x^-i to x^. The 
measurement made by the observer is an m-vector Zk, and is an m x n matrix that trans- 
forms a vector from the state space into the appropriate vector in the measurement space. The 
noise terms Wfe^y^-i (an n-vector) and (an m-vector) encompass known and unknown errors 
in Fk.k-\ and Hk and are normally distributed with mean and co variances given hy n x n 
matrix Qk,k-i and mxm matrix Rk, respectively. At each iteration, the Kalman Filter makes 
a state prediction for Xk, denoted Xj^\k-\- We use the notation k\k— 1 since we will only use 
measurements provided until time-step k—l in order to make the prediction at time-step k. 
The state prediction error Xk\k-\ is defined as the difference between the true state and the 
state prediction, as below. 

h\k-\=^k-h\k-\ (2.3) 

The covariance structure for the expected error on the state prediction is defined as the 
expectation of the outer product of the state prediction error. We call this covariance structure 
the error covariance prediction and denote it P^i^-i El 



k\k-l 



E 



{h\k-l) {h\k-l) 



(2.4) 



The filter will also provide an updated state estimate for Xk, given all the measurements 
provided up to and including time step k. We denote these estimates by Xk\k- We similarly 
define the state estimate error as below. 

^k\k = ^k-^k\k (2-5) 

The expectation of the outer product of the state estimate error represents the covariance 
structure of the expected errors on the state estimate, which we call the updated error covari- 
ance and denote P^i^. 



{h\k) {h\ky 



(2.6) 



At time-step k, we can make a prediction for the underlying state of the system by allowing 
the state to transition forward using our model for the dynamics and noting that E [wit^yt-i] = 0- 
This serves as our state prediction. 

h\k-\ = Fk,k-\Xk-\\k-\ (2-V) 

If we expand the expectation in Equation (|2.4I) . we have the following equation for the 
error covariance prediction. 

Pk\k-\ = Fk,k-\Pk-\\k-\Fk,k-\ + Qk,k-\ (2.8) 

We can transform our state prediction into the measurement space, which is a prediction 
for the measurement we now expect to observe. 



^We use the prime notation on a vector or a matrix to denote its transpose throughout this paper 
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h\k-i = HkXk\k-\ (2.9) 

The difference between the observed measurement and our predicted measurement is the 
measurement residual, which we are hoping to minimize in this algorithm. 

Vyt = Zfc-4|yt-l (2.10) 

We can also calculate the associated covariance for the measurement residual, which is the 
expectation of the outer product of the measurement residual with itself, E [v^v^] . We call this 
the measurement residual covariance. 

Sk = HkPk\k-iK + Rk (2.11) 

We can now define our updated state estimate as our prediction plus some perturbation, 
which is given by a weighting factor times the measurement residual. The weighting factor, 
called the Kalman Gain, will be discussed below. 

Xk\k=h\k-\+KkVk (2.12) 

Naturally, we can also calculate the updated error covariance by expanding the outer prod- 
uct in Equation (I2.6I) FI 

Pk\k = {l-KkHk)Pk\k-i {l-KM' + K^RkK (2.13) 

Now we would like to find the Kalman Gain Kj^, which minimizes the mean square state 
estimate error, E |%|^| . This is the same as minimizing the trace of the updated error 

covariance matrix aboveH After some calculus, we find the optimal gain that achieves this, 
written below@ 

Kk = Pk\k-iH',S^' (2.14) 

The covariance matrices in the Kalman Filter provide us with a measure for uncertainty 
in our predictions and updated state estimate. This is a very important feature for the various 
applications of filtering since we then know how much to trust our predictions and estimates. 
Also, since the method is recursive, we need to provide an initial covariance that is large 
enough to contain the initial state to ensure comprehensible performance. For a more detailed 
discussion of Kalman Filtering, we refer the reader to the following book [1]. 



"^The I in Equation (12.13b represents the nxn identity matrix. Throughout this paper, we use I to denote the 
same matrix, except in AppendixlAl where I is the appropriately sized identity matrix. 
^Note that v'v = trace [vv'] for some vector v. 

^We could also minimize the mean square state estimate error in the A' norm, where is a positive definite 
and symmetric weighting matrix. In the norm, the optimal gain would heK^ =N^Kk. 
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3 Equality Constrained Kalman Filtering 



A number of approaches have been proposed for solving the equality constrained Kalman 
Filtering problem [5, 11, 13, 14, 16]. In this paper, we show two different methods. The first 
method will restrict the state at each iteration to lie in the equality constrained space. The 
second method will start with a constrained prediction, and restrict the Kalman Gain so that 
the estimate will lie in the constrained space. Our equality constraints in this paper will be 
defined as below, where Ais a.qxn matrix, b a ^-vector, and x^, the state, is a n-vector|Zl 

Axk = b (3.1) 

So we would like our updated state estimate to satisfy the constraint at each iteration, as 
below. 

Axk\k = b (3.2) 

Similarly, we may also like the state prediction to be constrained, which would allow a 
better forecast for the system. 

Axk\k-i=b (3.3) 

In the following subsections, we will discuss methods for constraining the updated state 
estimate. In Section HI we will extend these concepts and formulations to the inequality con- 
strained case, and in Section [6l we will address the problem of constraining the prediction, as 
well. 

3.1 Projecting the state to lie in the constrained space 

We can solve the following minimization problem for a given time- step k, where is the 
constrained estimate, is any positive definite symmetric weighting matrix, and Xj^\j^ is the 
unconstrained Kalman Filter updated estimate. 

K\k = min \{x- 'Wk{x- -.Ax^b] (3.4) 
The best constrained estimate is then given by 

^k\k = ^k\k-W,-'A' {AW.'A'y' {Axk\k-b) (3.5) 

To find the updated error covariance matrix of the equality constrained filter, we first define 
the matrix T belowJl 

Y = Wj^^A' {AWj^^A') (3.6) 

Equation (13.51 ) can then be re- written as following. 

and b can be different for different k. We don't subscript each A and b to avoid confusion. 
**Note that TA is a projection matrix, as is (I — TA), by definition. If A is poorly conditioned, we can use a QR 
factorization to avoid squaring the condition number 
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We can find a reduced form for Xk — x^,,as below. 



(3.7) 



^k-xj^\k = ^k-Xk\k + '^ {Axk\k -b- {Axk - b)) (3.8a) 
= Xk-Xk\k + '^ {^h\k-A^k) (3.8b) 
= -{l-TA){xk\k-Xk) (3.8c) 

Using tlie definition of tlie error covariance matrix, we arrive at the following expression. 



p!\k=^ 


{^k-K\k) {^k-K\k) 




(3.9a) 


= E 


(I-TA) {xk\k-Xk) {xk\k 


-x,)'(I-TA)' 


(3.9b) 


= (I- 


-TA)P,|,(I-TA)' 




(3.9c) 


= Pk\k - TAPk\k - Pk\kAr + TAPk\kAr 


(3.9d) 


= Pj^^,-JAPkik 




(3.9e) 


= (I- 


-TA)P,|, 




(3.9f) 



It can be shown that choosing Wk = Pj^^^ results in the smallest updated error covariance. 
This also provides a measure of the information in the state at 



3.2 Restricting the optimal Kalman Gain so the updated state estimate 
lies in the constrained space 

Alternatively, we can expand the updated state estimate term in Equation (13.21) using Equation 
(EH. 

A{xkik-i+KkVk)=b (3.10) 

Then, we can choose a Kalman Gain K^, that forces the updated state estimate to be in 
the constrained space. In the unconstrained case, we chose the optimal Kalman Gain Kk, by 
solving the minimization problem below which yields Equation (12.141) . 

Kk = argmintrace [{l-KHk)Pk\k-i il-KHk)' + KRkK'] (3.11) 

^If M and are covariance matrices, we say is smaller than M if M — is positive semidefinite. Another 
formulation for incorporating equality constraints into a Kalman Filter is by observing the constraints as pseudo- 
measurements [14, 16]. When Wk is chosen to be P^l, both of these methods are mathematically equivalent [5]. 
Also, a more numerically stable form of Equation ( |3.9t with discussion is provided in [5]. 
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Now we seek the optimal that satisfies the constrained optimization problem written 
below for a given time-step k. 



= argmintrace [{\-KH,)P,,\k_, (l-KHk)' + KR,K'] 

KeRnxm (3.12) 

s.tA{xk\k-i+KVk)=b 

We will solve this problem using the method of Lagrange Multipliers. First, we take the 
steps below, using the vec notation (column stacking matrices so they appear as long vectors, 
see Appendix lA)) to convert all appearances of K in Equation (14.81) into long vectors. Let us 
begin by expanding the following term0 



trace [{l-KHk)Pk\k-i (l-KHk)' + KRkK'] 

= trace [Pk\j,_, -KHkP,\k , - Pk\k ,Kk' + KHkPk\k-iH[K' + KR^K'] 

™trace [Pk^^-i -KHj,P,\,_,-Pk\k-iH',K' + KSkK'] (3.13a) 
= trace [Pk\k-i] - trace [KHkPk\k^i] - trace [Pi,\k-iH'j^K'] + trace [KS^K'] 

We now expand the last three terms in Equation (|3.13a|) one at a timeHH 
trace [if/Z^P^in] ^ vec ' vec [K] 



(3.14) 



^wec[Pk^,_,H',] vec [if] 
trace [Pk\k-iH'kK'] ^ vec [if]'vec [Pk\k-iH'k\ (3.15) 



trace [KStK'] = vec [if]' vec [if5, 



(3.16) 



vec [if]' (5® I) vec [if] 

Remembering that trace [i\|jt-i] is constant, our objective function can be written as be- 
low. 

vec [if]' (I vec [if'] -vec [P^i.^i/i^j'vec [if] 

-vec[if]'vec[/',|,_i//^] 

Using Equation (IA.8I) on the equality constraints, our minimization problem is the follow- 
ing. 



'"Throughout this paper, a number in parentheses above an equals sign means we made use of this equation 
number. 

"We use the symmetry of Pm-i in Equation ( 13.141 1 and the symmetry of S/^ in Equation ( I3.16I I. 
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= arg min vec [K^ {Sk ® I) vec [K] 

-vec[/\|^_i//;]'vec[if] (3_jg) 

-vec[i^]'vec[/\|,_i//;] 

s.t. (v^ (g)A) vec [i^] = b -Ax^k-i 

Further, we simplify this problem so the minimization problem has only one quadratic 
term. We complete the square as follows. We want to find the unknown variable /i which will 
cancel the linear term. Let the quadratic term appear as follows. Note that the non-"vec [K]" 
term is dropped as is is irrelevant for the minimization problem. 

(vec [K]+^y {Sk ® I) (vec [K]+^) (3.19) 
The linear term in the expansion above is the following. 

vec [K]' {Sk ® I) M + Ai' (5^ ® I) vec [K] (3.20) 
So we require that the two equations below hold. 



(5^®I)M = -vec[%_i//;] 
li' {Sk0l) = -vec [Pk\k-iH'k]' 
This leads to the following value for /i. 



(3.21) 



Ai = -(5,i®l)vec[/\jfc_i//^] 
^ -vec [Pk\k-iH'kS,'] (3-22) 
= -vec[is:^] 

Using Equation (|A.6I) . our quadratic term in the minimization problem becomes the fol- 
lowing. 

(vec [K -Kk])' {Sk ® I) (vec [K - Kk] ) (3.23) 
Let / = vec [K — Kk] . Then our minimization problem becomes the following. 

is:^ = argmin l' {Sk0l)l 

leR'"" (3.24) 
s.t. {v'k(E)A){l + wec[Kk])=b-Axk\k-i 
We can then re-write the constraint taking the vec [Kk] term to the other side as below. 

(v^ ® A) l = b-Axk\k-i- {y'k ® A) vec [Kk] 



^ b -Axk\k-i - vec [AKkVk] 
= b-Axk\k-i-AKkVk 

= b-Axk\k 
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This results in the following simplified form. 

is:^ = argmin l' {Sk^l)l 

s.t. (v;®A)/ = Z>-A% 



(3.26) 



We form the Lagrangian where we introduce q Lagrange Multipliers in vector A 
(Ai, A2, • • ., A^) 



^=l'{Skm)l-?i'[ {v'k ®A)l-b +Mk\k] 
We take the partial derivative with respect to zl^ 



dl 



21' {Sk®l)-X' {\[®A) 



Similarly we can take the partial derivative with respect to the vector A. 



k\k 



(3.27) 



(3.28) 



(3.29) 



When both of these derivatives are set equal to the appropriate size zero vector, we have 
the solution to the system. Taking the transpose of Equation (13.281) . we can write this system 
as Mn = p with the following block definitions for M, n, and p. 



M 



25^® I Vk®A'' 







[qxq] 



^[mnx 1 

b-Ax 



(3.30) 
(3.31) 



(3.32) 



We solve this system for vector n in Appendix O The solution for I is pasted below. 



-1 



A'{AA') ' ){b-Axk\k) 



(3.33) 



Bearing in mind that b—Axj^^i^ = vec [Z? — Af^|^] , we can use Equation (IA.8I) to re- write / 
as belowll^ 



vec 



A' [AA')-' {b-Ax,\,) {y'kS,'vky' y'kS,' 



-1 



(3.34) 



The resulting matrix inside the vec operation is then an n by m matrix. Remembering the 
definition for /, we notice that K — Kj^ results in an n by m matrix also. Since both of the 
components inside the vec operation result in matrices of the same size, we can safely remove 



12 



13 



We used the symmetry of (8^ 1) here. 

Here we used the symmetry of 5^' and {vlS^^Vk) ' (the latter of which is actually just a scalar). 
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the vec operation from both sides. This results in the following optimal constrained Kalman 
Gain K^. 

K,,-A'{AA!y' {Axk\^-b) {v',S,'vky\',Sj;' (3.35) 

If we now substitute this Kalman Gain into Equation (12.121) to find the constrained updated 
state estimate, we end up with the following. 

^k\k=^k\k-A'{AA'y' {Ax^k-b) (3.36) 

This is of course equivalent to the result of Equation (|3.5I) with the weighting matrix Wk 
chosen as the identity matrix. The error covariance for this estimate is given by Equation 



4 Adding Inequality Constraints 

In the more general case of this problem, we may encounter equality and inequality constraints, 
as given belowl^ 

Axk = b 

(4.1) 

Cxic < d 

So we would like our updated state estimate to satisfy the constraint at each iteration, as 
below. 

Axui. = b 

. (4.2) 

Cxj^^j^ ^ d 

Similarly, we may also like the state prediction to be constrained, which would allow a 
better forecast for the system. 

Axui,_, = b 

k\k ^^^^ 

Ciyt|yfc-1 < d 

We will present two analogous methods to those presented for the equality constrained 
case. In the first method, we will run the unconstrained filter, and at each iteration constrain 
the updated state estimate to lie in the constrained space. In the second method, we will 
find a Kalman Gain such that the the updated state estimate will be forced to lie in the 
constrained space. In both methods, we will no longer be able to find an analytic solution as 
before. Instead, we use numerical methods. 



'^We can use the unconstrained or constrained Kalman Gain to find this error covariance matrix. Since the 
constrained Kalman Gain is suboptimal for the unconstrained problem, before projecting onto the constrained 
space, the constrained covariance will be different from the unconstrained covariance. However, the difference 
lies exactly in the space orthogonal to which the covariance is projected onto by Equation (13.9b . The proof is 
omitted for brevity. 

'^C and d can be different for different k. We don't subscript each C and d to simplify notation. 
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4.1 By Projecting the Unconstrained Estimate 



Given the best unconstrained estimate, we could solve the following minimization problem 
for a given time-step k, where is the inequality constrained estimate and is any positive 
definite symmetric weighting matrix. 

•^A: = argmin {x - x^k)' Wk {x - Xk^k) 

It. Ax = b (4.4) 
Cx < d 

For solving this inequality constrained optimization problem, we can use a variety of stan- 
dard methods, or even an out-of-the-box solver, like f mincon in Matlab. Here we use an ac- 
tive set method [4]. This is a common method for dealing with inequality constraints, where 
we treat a subset of the constraints (called the active set) as additional equality constraints. 
We ignore any inactive constraints when solving our optimization problem. After solving the 
problem, we check if our solution lies in the space given by the inequality constraints. If it 
doesn't, we start from the solution in our previous iteration and move in the direction of the 
new solution until we hit a set of constraints. For each iteration, the active set is made up of 
those inequality constraints with non-zero Lagrange Multipliers. 

We first find the best estimate (using Equation (|3.5I) for the equality constrained problem 
with the equality constraints given in Equation (14.11 ) plus the active set of inequality constraints. 
Let us call the solution to this xf^ ■ since we have not yet checked if the solution lies in the 

inequality constrained spacel^ In order to check this, we find the vector that we moved along 
to reach i^j^ j. This is given by the following. 

^=^\k,j-K\k,j~i (4-5) 

We now iterate through each of our inequality constraints, to check if they are satisfied. If 
they are all satisfied, we choose Tmax = 1- If they are not, we choose the largest value of Tmax 
such that iyt|A;j-i + '^max.s lics in the inequality constrained space. We choose our estimate to 
be 



^\k,j ~ ^\kj-l + '^max* (4.6) 

If we find the solution has converged within a pre-specified error, or we have reached a 
pre-specified maximum number of iterations, we choose this as the updated state estimate to 
our inequality constrained problem, denoted i^^. If we would like to take a further iteration 

on j, we check the Lagrange Multipliers at this new solution to determine the new active set0 
We then repeat by finding the best estimate for the equality constrained problem including 
the new active set as additional equality constraints. Since this is a Quadratic Programming 
problem, each step of j guarantees the same estimate or a better estimate. 

'^For the inequality constrained filter, we allow multiple iterations within each step. The j subscript indexes 
these further iterations. 

'^The previous active set is not relevant. 
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When calculating the error covariance matrix for this estimate, we can also add on the 
safety term below. 



This is a measure of our convergence error and should typically be small relative to the 
unconstrained error covariance. We can then use Equation (13.91) to project the covariance 
matrix onto the constrained subspace, but we only use the defined equality constraints. We do 
not incorporate any constraints in the active set when computing Equation (13.91) since these 
still represent inequality constraints on the state. Ideally we would project the error covariance 
matrix into the inequality constrained subspace, but this projection is not trivial. 

4.2 By Restricting the Optimal Kalman Gain 

We could solve this problem by restricting the optimal Kalman gain also, as we did for equality 
constraints previously. We seek the optimal K]^ that satisfies the constrained optimization 
problem written below for a given time- step k. 



Again, we can solve this problem using any inequality constrained optimization method 
(e.g., fmincon in Matlab or the active set method used previously). Here we solved the 
optimization problem using SDPT3, a Matlab package for solving semidefinite programming 
problems [15]. When calculating the covariance matrix for the inequality constrained estimate, 
we use the restricted Kalman Gain. Again, we can add on the safety term for the convergence 
error, by taking the outer product of the difference between the updated state estimates cal- 
culated by the restricted Kalman Gain for the last two iterations of SDPT3. This covariance 
matrix is then projected onto the subspace as in Equation (13.91 ) using the equality constraints 
only. 

5 Dealing with Nonlinearities 

Thus far, in the Kalman Filter we have dealt with linear models and constraints. A number 
of methods have been proposed to handle nonlinear models (e.g.. Extended Kalman Filter [1], 
Unscented Kalman Filter [7]). In this paper, we will focus on the most widely used of these, the 
Extended Kalman Filter. Let's re-write the discrete unconstrained Kalman Filtering problem 
from Equations (12. 1|) and (12.21) below, incorporating nonlinear models. 




(4.7) 



kl = argmintrace [{l-KHk)Pk\k-i (l-KHk)' + KRkK'] 



s.tA{xk\k_i+KkVk) =b 
C {k\k-i + Kkyk) <d 



(4.8) 



— fk,k-\ {^k-l) +l^k,k-li 



Uk,k-\ ~ ^ (0,2yt,yt-l) 



(5.1) 



Zk = hk {xk)+Vk., 



(5.2) 
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In the above equations, we see that the transition matrix i has been replaced by the 
nonlinear vector-valued function /^ yt-i (■)' similarly, the matrix Hj^, which transforms a 
vector from the state space into the measurement space, has been replaced by the nonlinear 
vector- valued function hk{-). The method proposed by the Extended Kalman Filter is to lin- 
earize the nonlinearities about the current state prediction (or estimate). That is, we choose 
Fk,k-\ as the Jacobian of fk,k-i evaluated at i;t-i|yt-i ' ^k as the Jacobian of % evaluated at 
and proceed as in the linear Kalman Filter of Section [20 Numerical accuracy of these 
methods tends to depend heavily on the nonlinear functions. If we have linear constraints 
but a nonlinear fk.k-i (•) and hk{-), we can adapt the Extended Kalman Filter to fit into the 
framework of the methods described thus far. 



5.1 Nonlinear Equality and Inequality Constraints 

Since equality and inequality constraints we model are often times nonlinear, it is important to 
make the extension to nonlinear equality and inequality constrained Kalman Filtering for the 
methods discussed thus far. Without loss of generality, our discussion here will pertain only to 
nonlinear inequality constraints. We can follow the same steps for equality constraints We 
replace the linear inequality constraint on the state space by the following nonlinear inequal- 
ity constraint c{xk) = d, where c(-) is a vector-valued function. We can then linearize our 
constraint, c (xk) = d, about the current state prediction Xj^\k-i^ which gives us the following^ 

c{xk\k-i)+C{xk-Xkik_i) ^d (5.3) 

Here C is defined as the Jacobian of c evaluated at ifc|A:-i- This indicates then, that the 
nonlinear constraint we would like to model can be approximated by the following linear 
constraint 



Cxk^d + Cxi,\,,_i-c{xk\k_i) (5.4) 

This constraint can be written as Cx^ < d, which is an approximation to the nonlinear 
inequality constraint. It is now in a form that can be used by the methods described thus far. 

The nonlinearities in both the constraints and the models, fk,k-\ (■) and /jyt (■), could have 
been linearized using a number of different methods (e.g., a derivative-free method, a higher 
order Taylor approximation). Also an iterative method could be used as in the Iterated Ex- 
tended Kalman Filter [1]. 



^'^We can also do a midpoint approximation to find Fi^j^^x by evaluating the Jacobian at +x^i^_\) /2. 

This should be a much closer approximation to the nonlinear function. We use this approximation for the Ex- 
tended Kalman Filter experiments later. 

'^We replace the '<' sign with an sign and the with an sign. 

^"This method is how the Extended Kalman Filter linearizes nonlinear functions for fk,k-\ (•) and /i/i (•). Here 
ijtlit-i can be the state prediction of any of the constrained filters presented thus far and does not necessarily relate 
to the unconstrained state prediction. 
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6 Constraining the State Prediction 



We haven't yet discussed whether the state prediction (Equation (12.71) ) also should be con- 
strained. Forcing the constraints should provide a better prediction (which is used for fore- 
casting in the Kalman Filter). Ideally, the transition matrix F^ ^^-i will take an updated state 
estimate satisfying the constraints at time k—\ and make a prediction that will satisfy the con- 
straints at time k. Of course this may not be the case. In fact, the constraints may depend on the 
updated state estimate, which would be the case for nonlinear constraints. On the downside, 
constraining the state prediction increases computational cost per iteration. 

We propose three methods for dealing with the problem of constraining the state predic- 
tion. The first method is to project the matrix F^^k-i onto the constrained space. This is only 
possible for the equality constraints, as there is no trivial way to project F^^k-i to an inequal- 
ity constrained space. We can use the same projector as in Equation (|3.9f|) so we have the 
following 



^L^-i = (I-Tf^)^U-i (6.1) 

Under the assumption that we have constrained our updated state estimate, this new transi- 
tion matrix will make a prediction that will keep the estimate in the equality constrained space. 
Alternatively, if we weaken this assumption, i.e., we are not constraining the updated state 
estimate, we could solve the minimization problem below (analogous to Equation (13.41 )). We 
can also incorporate inequality constraints now. 

^k\k- 1 = arg min {x - i ) ' Wj, {x - Xk\k- 1 ) 

X 

s.t. Ax = b (6.2) 

Cx <d 

We can constrain the covariance matrix here also, in a similar fashion to the method de- 
scribed in Section 14.11 The third method is to add to the constrained problem the additional 
constraints below, which ensure that the chosen estimate will produce a prediction at the next 
iteration that is also constrained. 

^k+iFk+i^k^k = bk+1 ^ , „^ 

{o.j) 

Ck+iFk+i,k-Xk < dk+i 

If Ayt+i^^/t+iiQ+i or depend on the estimate (e.g., if we are linearizing nonlinear 
functions a(-) or b{-), we can use an iterative method, which would resolve A^t+i and b^^i 
using the current best updated state estimate (or prediction), re-calculate the best estimate 
using Aj^^i and b^^i, and so forth until we are satisfied with the convergence. This method 
would be preferred since it looks ahead one time-step to choose a better estimate for the current 
iteration j^^l However, it can be far more expensive computationally. 

■^'in these three methods, the symmetric weighting matrix can be different. The resulting T can conse- 
quently also be different. 

^^Further, we can add constraints for some arbitrary n time-steps ahead. 
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7 Experiments 



We provide two related experiments here. We have a car driving along a straight road with 
thickness 2 meters. The driver of the car traces a noisy sine curve (with the noise lying only in 
the frequency domain). The car is tagged with a device that transmits the location within some 
known error. We would like to track the position of the car. In the first experiment, we filter 
over the noisy data with the knowledge that the underlying function is a noisy sine curve. The 
inequality constrained methods will constrain the estimates to only take values in the interval 
[—1,1]. In the second experiment, we do not use the knowledge that the underlying curve is a 
sine curve. Instead we attempt to recover the true data using an autoregressive model of order 
6 [3]. We do, however, assume our unknown function only takes values in the interval [—1,1], 
and we can again enforce these constraints when using the inequality constrained filter. 

The driver's path is generated using the nonlinear stochastic process given by Equation 
(|5.1I) . We start with the following initial point. 



Xo 



Om 
Om 



(7.1) 



Our vector- valued transition function will depend on a discretization parameter T and can 
be expressed as below. Here, we choose T to be n/lO, and we run the experiment from an 
initial time of to a final time of lOn. 



k,k-l 



{Xk-l)i+T 

sm{{xk-i)i + T) 



(7.2) 



And for the process noise we choose the following. 



Qk,k-\ 



0.1 m^ 
0m2 



(7.3) 



The driver's path is drawn out by the second element of the vector x^ - the first element 
acts as an underlying state to generate the second element, which also allows a natural method 
to add noise in the frequency domain of the sine curve while keeping the process recursively 
generated. 



7.1 First Experiment 



To create the measurements, we use the model from Equation (12.21 ), where is the square 
identity matrix of dimension 2. We choose Rk as below to noise the data. This considerably 
masks the true underlying data as can be seen in Fig. fTj^ 



Rk 



10 m2 






10 m^ 



(7.4) 



^^The figure only shows the noisy sine curve, which is the second element of the measurement vector. The 
first element, which is a noisy straight line, isn't plotted. 
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Figure 1: We take our sine curve, which is already noisy in the frequency domain (due to process noise), 
and add measurement noise. The underlying sine curve is significantly masked. 



For the initial point of our filters, we choose the following point, which is different from 
the true initial point given in Equation (17.11) . 



Om 
1 m 



(7.5) 



Our initial covariance is given as below@. 

_rim2 0.1 
0.1 lm2 



(7.6) 



^^Nonzero off-diagonal elements in the initial covariance matrix often help the filter converge more quickly 
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In the filtering, we use the information that the underlying function is a sine curve, and our 
transition function fk^k-\ changes to reflect a recursion in the second element of - now we 
will add on discretized pieces of a sine curve to our previous estimate. The function is given 
explicitly below. 



k,k-l 



{xk-i)i+sm{{xk-i)i+T)-sm{{xk-i)i] 



(7.7) 



For the Extended Kalman Filter formulation, we will also require the Jacobian of this 
matrix denoted Fk_k-\, which is given below. 



k±-\ 



1 



cos((x^_i)i + r)-cos((x^_i)i; 



(7.8) 



The process noise Qk,k-\^ given below, is chosen similar to the noise used in generating 
the simulation, but is slightly larger to encompass both the noise in our above model and to 
prevent divergence due to numerical roundoff errors. The measurement noise Rk is chosen the 
same as in Equation (|7.4I) . 



Qk,k- 



0.1 m2 
0.1 m2 



(7.9) 



The inequality constraints we enforce can be expressed using the notation throughout the 
chapter, with C and d as given below. 



C 



1 
-1 



(7.10) 



(7.11) 



These constraints force the second element of the estimate (the sine portion) to lie in 
the interval [—1,1]. We do not have any equality constraints in this experiment. We run the 
unconstrained Kalman Filter and both of the constrained methods discussed previously. A plot 
of the true position and estimates is given in Fig.[2l Notice that both constrained methods force 
the estimate to lie within the constrained space, while the unconstrained method can violate 
the constraints. 



7.2 Second Experiment 

In the previous experiment, we used the knowledge that the underlying function was a noisy 
sine curve. If this is not known, we face a significantly harder estimation problem. Let us 
assume nothing about the underlying function except that it must take values in the interval 
[—1,1]. A good model for estimating such an unknown function could be an autoregressive 
model. We can compare the unconstrained filter to the two constrained methods again using 
these assumption and an autoregressive model of order 6, or AR(6) as it is more commonly 
referred to. 
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True Position and Estimates 




5 10 15 20 25 30 

Time (seconds) 



Figure 2: We show our true underlying state, which is a sine curve noised in the frequency domain, 
along with the estimates from the unconstrained Kalman Filter, and both of our inequality constrained 
modifications. We also plotted dotted horizontal lines at the values -1 and 1. Both inequality con- 
strained methods do not allow the estimate to leave the constrained space. 

In the previous example, we used a large measurement noise Rk to emphasize the gain 
achieved by using the constraint information. Such a large is probably not very realistic, 
and when using an autoregressive model, it will be hard to track such a noisy signal. To 
generate the measurements, we again use Equation (12.21) . this time with //^ and R^ as given 
below. 

Hk=[Q 1] (7.12) 
Rk=[0.5m^] (7.13) 
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Our state will now be defined using the following 13-vector, in which the first element is 
the current estimate, the next five elements are lags, the six elements afterwards are coefficients 
on the current estimate and the lags, and the last element is a constant term. 

^k\k=[yk Jk-i ■■■ yk-5 «! CC2 ■■■ a-j]' (7.14) 

Our matrix in the filter is a row vector with the first element 1, and all the rest as 0, so 
y^^k-i is actually our prediction Zk\k-i i^i the filter, describing where we believe the expected 
value of the next point in the time-series to lie. For the initial state, we choose a vector of all 
zeros, except the first and seventh element, which we choose as 1 . This choice for the initial 
conditions leads to the first prediction on the time series being 1 , which is incorrect as the true 
underlying state has expectation 0. For the initial covariance, we choose I[i3xi3] and add 0.1 
to all the off-diagonal elements El The transition function fk^k-i for the AR(6) model is given 
below. 

min(l,max(-l,ai3;;t-i H h tteJ/t-e + «?)) 

min(l,max(-l,3;^_i)) 

min(l,max(-l,3;^_2)) 
min(l,max(-l,3;^_3)) 
min(l,max(-l,j^_4)) 
min(l,max(-l,y^_5)) 
oci 

"2 



06 

Putting this into recursive notation, we have the following. 

+ (■^^-1)13)) 



(7.15) 



min(l,max(— 1, (jCyt-i)? 


{Xk-l)^ + 


...4 


min ( 1 , max ( 


-^:{Xk-l 


)l)) 


min ( 1 , max ( 


-l,(^yt-l 


)2)) 


min ( 1 , max ( 


-l,(-«yt-l 


)3)) 


min ( 1 , max ( 


-i^iXk-l 


)4)) 


min ( 1 , max ( 


-i^ixk-i 


)5)) 


{Xk- 


-1)7 




(Xk- 


-1)8 




[Xk- 


1)12 




{Xk- 


1)13 





(7.16) 



The Jacobian of fk^k-i is given below. We ignore the min(-) and max(-) operators since 
the derivative is not continuous across them, and we can reach the bounds by numerical error. 



^^The bracket subscript notation is used through the remainder of this paper to indicate the size of zero matrices 
and identity matrices. 
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Further, when enforced, the derivative would be 0, so by ignoring them, we are allowing our 
covariance matrix to be larger than necessary as well as more numerically stable. 



[5x5] 







[5x1] 



Oi 



7x6] 



0[5x7] 
i[7x7] 



(7.17) 



For the process noise, we choose Qk,k-i to be a diagonal matrix with the first entry as 0.1 
and all remaining entries as 10^^ since we know the prediction phase of the autoregressive 
model very well. The inequality constraints we enforce can be expressed using the notation 
throughout the chapter, with C as given below and J as a 12- vector of ones. 



C 



I 



[6x6] 



l[6x6] 



Oi 



12x7] 



(7.18) 



These constraints force the current estimate and all of the lags to take values in the range 
[— 1 , 1] . As an added feature of this filter, we are also estimating the lags at each iteration using 
more information although we don't use it - this is a fixed interval smoothing. In Fig. [3l we 
plot the noisy measurements, true underlying state, and the filter estimates. Notice again that 
the constrained methods keep the estimates in the constrained space. Visually, we can see the 
improvement particularly near the edges of the constrained space. 



8 Conclusions 

We've provided two different formulations for including constraints into a Kalman Filter. In 
the equality constrained framework, these formulations have analytic formulas, one of which 
is a special case of the other. In the inequality constrained case, we've shown two numerical 
methods for constraining the estimate. We also discussed how to constrain the state prediction 
and how to handle nonlinearities. Our two examples show that these methods ensure the 
estimate lies in the constrained space, which provides a better estimate structure. 



Appendix A Kron and Vec 



In this appendix, we provide some definitions used earlier in the chapter. Given matrix A G 
j^raxn g ^ M^^"?, we can define the right Kronecker product as belowj^ 



(A® 5) 



ai.i5 



a\nB 



^m,nB 



(A.l) 



^^The indices m,n,p, and q and all matrix definitions are independent of any used earlier. Also, the subscript 
notation ai„ denotes the element in the first row and n-th column of A, and so forth. 
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Measurement, True Position, and Estimates 




5 10 15 20 25 30 



Time (seconds) 

Figure 3: We show our true underlying state, which is a sine curve noised in the frequency domain, the 
noised measurements, and the estimates from the unconstrained and both inequality constrained filters. 
We also plotted dotted horizontal lines at the values -1 and 1. Both inequality constrained methods do 
not allow the estimate to leave the constrained space. 

Given appropriately sized matrices A,5,C, and D such that all operations below are well- 
defined, we have the following equalities. 

{A®B)' ={A'®B') (A.2) 
{A®By^ = [A-^®B-^) (A.3) 
(A ® 5) (C ® D) = (AC ® BD) (A.4) 
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We can also define the vectorization of an [m x n] matrix A, wliicli is a linear transformation 
on a matrix that stacks the columns iteratively to form a long vector of size \mn x 1], as below. 



«1,2 

ai,n 



vec [A] 



Using the vec operator, we can state the trivial definition below. 

vec [A + 5] = vec [A] + vec [B] 
Combining the vec operator with the Kronecker product, we have the following. 

vec[A5] = (5'® I) vec [A] 



(A.5) 



(A.6) 



(A.7) 



vec [ABC] = {C' ®A) vec [B] 
We can express the trace of a product of matrices as below. 



(A.8) 



trace [AB] = vec [B'] ' vec [A] 



(A.9) 



trace [ABC] = vec [5]' (I ®C) vec [A] (A. 10a) 

= vec[A]'(I®5)vec[C] (A.lOb) 

= vec[A]'(C®I)vec[5] (A. 10c) 

For more information, please see [9]. 



Appendix B Analytic Block Representation for the inverse 
of a Saddle Point Matrix 

Ms is a saddle point matrix if it has the block form belowEZl 

^^The subscript S notation is used to differentiate these matrices from any matrices defined earlier 
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Ms 



As B's 
Bs -Cs 



(B.l) 



In the case that As is nonsingular and the Schur complement Js = — {'^s + BsAg^B'^) is 
also nonsingular in the above equation, it is known that the inverse of this saddle point matrix 
can be expressed analytically by the following equation (see e.g., [2]). 



M 



-1 



Ag 1 +Ag ^B'gJg ^BsAg 1 -Ag ^B'gJ^ ^ 
1d_a 1 j-i 



Appendix C Solution to the system Mn = p 



(B.2) 



Here we solve the system Mn = p from Equations (|3.30l) . (13.311) . and (13.321) . re-stated below, 
for vector n. 



vi^A 



qxq] 

















p-AXklk_ 



(C.l) 



M is a saddle point matrix with the following equations to fit the block structure of Equa- 



tion dEB 



We can calculate the term A^ ^B'g. 



As = 2Sk®l 
Bs = Vk®A 

Cs = %xq] 



A-s'B^s=[USk®W' 



1a31a31i 1 



(5^-1® I) {vk®A') 



Ell 1 



(S,Vfc) ®A' 



And as a result we have the following for Js- 



(C.2) 
(C.3) 
(C.4) 



(C.5a) 
(C.5b) 

(C.5c) 



^-\Hs,\,)®{aa') 



(C.6a) 
(C.6b) 



28 



We use Equation ( IA.2b with B'^ to arrive at the same term for Bs in Equation ( IC.ll l. 
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is then, as below. 



7 



For the upper right block of M , we then have the following expression. 



(C.7a) 
(C.7b) 



-1 



Sk'^kiykS^Vk) ® A'{AA') 



(C.8a) 
(C.8b) 



Since the first block element of p is a vector of zeros, we can solve for n to arrive at the 
following solution for /. 



Sk ^kiVkSyVk) ® A'{AA') ){b-Axkik) 



The vector of Lagrange Multipliers A is given below. 



-2 



(C.9) 



(C.IO) 



References 

[1] Y. Bar-Shalom, X. R. Li, and T. Kirubarajan, Estimation with Applications to 
Tracking and Navigation, John Wiley and Sons, Inc., 2001. 

[2] M. Benzi, G. H. Golub, and J. Liesen, Numerical solution of saddle point problems. 
Acta Numerica, 14 (2005), pp. 1-137. 

[3] G. E. P. Box AND G. M. Jenkins, Time Series Analysis. Forecasting and Control 
(Revised Edition), Oakland: Holden-Day, (1976). 

[4] R. Fletcher, Practical methods of optimization. Vol. 2: Constrained Optimiza- 
tion, John Wiley & Sons Ltd., Chichester, 1981. Constrained optimization, A Wiley- 
Interscience Publication. 



[5] N. Gupta, Kalman filtering in the presence of state space equality constraints, in 
IEEE Proceedings of the 26th Chinese Control Conference, July 2007, arXiv:physics.ao- 
ph/0705.4563, Oxford na-tr:07/14. 

[6] N. Gupta, R. Hauser, and N. F. Johnson, Forecasting financial time-series using 
an artificial market model, in Proceedings of the 10th Annual Workshop on Economic 
Heterogeneous Interacting Agents, June 2005, Oxford na-tr:05/09. 



25 



[7] S. J. JULIER AND J. K. Uhlmann, A new extension of the kalman filter to non- 
linear systems, in Proceedings of AeroSense: The 11th International Symposium on 
Aerospace/Defence Sensing, Simulation and Controls, vol. 3, 1997, pp. 182-193. 

[8] R. E. Kalman, A new approach to linear filtering and prediction problems. Transac- 
tions of the ASME- Journal of Basic Engineering, 82 (1960), pp. 35^5. 

[9] P. Lancaster and M. Tismenetsky, The Theory of Matrices: With Applications, 
Academic Press Canada, 1985. 

[10] A. G. QURESHI, Constrained kalman filtering for image restoration, in Proceedings of 
the International Conference on Acoustics, Speech, and Signal Processing, vol. 3, 1989, 
pp. 1405 - 1408. 

[11] D. Simon and T. L. Chia, Kalman filtering with state equality constraints, IEEE 
Transactions on Aerospace and Electronic Systems, 38 (2002), pp. 128-136. 

[12] D. Simon and D. L. Simon, Aircraft turbofan engine health estimation using con- 
strained kalman filtering. Journal of Engineering for Gas Turbines and Power, 127 
(2005), p. 323. 

[13] T. L. Song, J. Y. Ahn, and C. Park, Suboptimal filter design with pseudomeasure- 
mentsfor target tracking, IEEE Transactions on Aerospace and Electronic Systems, 24 
(1988), pp. 28-39. 

[14] M. Tahk and J. L. Speyer, Target tracking problems subject to kinematic constraints. 
Proceedings of the 27th IEEE Conference on Decision and Control, (1988), pp. 1058- 
1059. 

[15] K. C. TOH, M. J. Todd, and R. Tutuncu, SDPT3 — a Matlab software package 
for semidefinite programming. Optimization Methods and Software, 11 (1999), pp. 545- 
581. 

[16] L. S. Wang, Y. T. Chiang, and F. R. Chang, Filtering method for nonlinear sys- 
tems with constraints, lEE Proceedings - Control Theory and Applications, 149 (2002), 
pp. 525-531. 



26 



